A Bell-Evans-Polanyi principle for molecular dynamics trajectories and its 

implications for global optimization 
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The Bell-Evans-Polanyi principle that is valid for a chemical reaction that proceeds along the 
reaction coordinate over the transition state is extended to molecular dynamics trajectories that 
in general do not cross the dividing surface between the initial and the final local minima at the 
exact transition state. Our molecular dynamics Bell-Evans-Polanyi principle states that low energy 
molecular dynamics trajectories are more likely to lead into the basin of attraction of a low energy 
local minimum than high energy trajectories. In the context of global optimization schemes based 
on molecular dynamics our molecular dynamics Bell-Evans-Polanyi principle implies that using low 
energy trajectories one needs to visit a smaller number of distinguishable local minima before finding 
the global minimum than when using high energy trajectories. 



The Bell-Evans-Polanyi (BEP) principle [J Q is an 
important fundamental principle in chemistry. It gives a 
relation between the free energy AG released in a chem- 
ical reaction and the activation free energy e a for the re- 
action. It was qualitatively first put forward by Br0nsted 
•3ij who observed that strongly exothermic reactions have 
a low activation energy. A more quantitative relation was 
then derived by Polanyi et a/[l], Q who approximated the 
potential energy surface by straight lines. This approxi- 
mation leads to a linear relation between the activation 
energy e a and the free energy of the reaction AG: 
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fci + fc 2 AG 



(1) 



where k\ and hi are constants that depend on the slopes 
of the lines. A more accurate approach by Marcus [J Q 
approximates the potential energy surface by two parabo- 
las centered at the two local minima of the energy. The 
potential energy surface in this approximation is every- 
where a quadratic form with the exception of the inter- 
section point of the two parabolas where it has a discon- 
tinuity in its derivative. From FigQ] it is easily visible 
that the barrier for the reaction A — > B is lowered if the 
parabola centered in the minimum B is shifted downward. 
The resulting quantitative relation^ is given by 




AG AG 2 
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(2) 



where k is proportional to the curvature of the two 
parabolas. 

In a chemical reaction the reaction coordinate con- 
nects the educt A with the product B. Hence the in- 
tersection of the two parabolas is the transition state. 
In this article we will study the BEP principle not for 
this hypothetical path along the reaction coordinate but 
for molecular dynamics (MD) trajectories that cross the 
dividing hypersurface between the two basins of attrac- 
tion Q of two local minima on the potential energy sur- 
face. The notions of educt and product are replaced by 
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FIG. 1: The potential energy surface in the region between 
minimum A and B in the Marcus approximation. The poten- 
tial energy in the neighborhood of minimum A is given by a 
parabola centered at A, and in the neighborhood of minimum 
B by a parabola centered at B. The activation energy e a for 
the chemical reaction A — > B is determined by the intersec- 
tion of the two parabolas. It can be seen that the activation 
energies are smaller if the local minimum B is lower in energy. 



the notions of initial and final local minima in this con- 
text. We will show that the BEP principle is also valid 
in the context of MD. 

Since our study requires the calculation and statistical 
evaluation of a very large number of local minima and 
saddle points, we will initially base our study on a 
Lennard Jones cluster [t| containing 55 atoms for 
which stationary points can be calculated rapidly. The 
parameters entering in the Lennard Jones potential were 
selected such that it models Argon clusters, namely 
e=0.998 kJ/Mole, a=3AA and M=39.948amu Q. The 
free energies were calculated within the harmonic ap- 
proximation as the vibrational free energy. The tern- 
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perature at which the free energies were calculated is 
30 K which is below the melting point (50 K) of this 
weakly bound system 0(Lennard Jones cluster). Ini- 
tially we have searched for more than 130000 first order 
saddle points Gf on the potential energy surface con- 
necting energetically low local minima. Subsequently we 
have moved the system to the left and to the right along 
the direction where the curvature is negative. These two 
points served as the starting points for a local geometry 
optimization that led us in the two closest local minima. 
In this way we have generated pairs of local minima to- 
gether with the saddle points that connect them. Fig. [3] 
and Fig. [2] show scatter plots of AG = G\ — Gf versus 
the activation energy e a — Gf — Gf with and without the 
entropy contributions respectively and Fig.[4]shows a his- 
togram with averages of the Gf — Gf. Each pair of local 
minima contributed two data points to these plots since 
one can surmount the barrier by going from the minimum 
A to minimum B as well as by going from minimum B 
to minimum A. 




Ac:iva:ion Energy (Epsibn unh) 

FIG. 2: The relation between the activation energy G s — G a 
and the reaction energy G\ — Gf for more than 130000 saddle 
points in a Lennard Jones cluster of 55 atoms. All the energies 
plotted here are free energies at T = 0, i.e. just energies 



The two scatter plots in Fig. [2] and Fig. [3] show that 
there is no strict linear correlation between the barrier 
height e a and the energy difference AG between the two 
minima. For small barrier heights one can find both high 
energy and low energy minima behind the barrier. How- 
ever, the BEP principle holds as a negation. If one goes 
over high barriers it is extremely unlikely that one will 
end up in a low energy minimum. The better correla- 
tion for large activation energies is simply due to the fact 
that AG can not become larger than e a . On the other 
hand, Fig. [4] shows that there is a good linear relation if 
one averages over AG. Good linear Bell-Evans-Polanyi 
relations have been found in calculations of dissociative 
chemisorption of various molecules [H Q M, [3 . 
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FIG. 3: The relation between the activation energy G s 
and the reaction energy G\ 
points at T = 30 K. 
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FIG. 4: The same data as in Fig. [2] and Fig. [3] but averaged 
within 25 bins along the x axis. 



Kinetic rate theory gives the rate constant for a reac- 
tion as 

k = !^exp{-0e a ) = ^^exp{P(E*-E>)) , 

where E a and Eb are the energies of the two minima. Q s 
is the partition function for the transition state and gives 
in a certain sense the size of the barrier region. The im- 
portance of entropy terms can easily be seen in the clas- 
sical limit. By making the same approximation as was 
done originally by Marcus in the derivation of the BEP 
principle, namely that the potential energy surface is 
the union of parabolas, but by considering 2-dimcnsional 
parabolas instead of 1-dimensional parabolas, one can 
easily see in Fig. [5] that the size of the crossing surface 
that can be surmounted by a MD path of limited energy 
is increasing strongly when the MD path goes into an 
energetically low basin. We expect therefore that for a 
fixed energy crossings into low energy minima are more 
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probable than crossings into high energy local minima. 
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FIG. 6: The number of visits as a function of Et — E a 



FIG. 5: This figure shows in blue the potential energy sur- 
face of 3 local minima together with their basins of attrac- 
tion. The value at the minimum is for the minimum in 
the middle and in the background, but -0.5 for the minimum 
in the foreground. The region where a MD trajectory of en- 
ergy .25 can cross from the central basin of attraction into 
the foremost basin of attraction (green) is much wider than 
the region where the trajectory can cross into the basin in 
the background (red). The value of the partition function Q s 
for crossings into the basin in the foreground is thus expected 
to be much larger than for crossings into the basin in the 
background. 



Fig. [S] shows the results of a numerical experiment. 
For MD trajectories that start with random directions 
but fixed kinetic energy Ekin from a certain minimum 
with energy E a we have recorded how many times this 
trajectory reaches the basin of attraction of neighbor- 
ing minima with energy E^. To check whether the MD 
trajectory has crossed into another basin of attraction 
steepest descent geometry optimizations were started af- 
ter every 20 MD steps. In Fig. [6] we then plot the number 
of visits as a function of Eb — E a . We see that it is orders 
of magnitude more likely that the MD trajectory crosses 
into low energy basins than in high energy basins. We 
will denote this correlation as the MDBEP principle: low 
energy MD trajectories are more likely to lead into the 
basin of attraction of a low energy local minimum than 
high energy trajectories. The activation energy of the 
original BEP principle has thus been replaced by the en- 
ergy of the trajectory. This implies that we have replaced 
of property of the potential energy surface by a property 
of the MD trajectory exploring this surface. 

As can be seen from Fig. [2 Fig. [3] and FiglH both the 
traditional BEP principle and our MDBEP principle are 
only valid in an average sense. Such a validity on the 
average is sufficient in the context of global optimiza- 
tion using the minima hopping method (MHM) [1, 0]. 
In MHM the system moves from one local minimum to 
another by a combination of MD followed by a local ge- 



ometry optimization. With the MD part one jumps from 
one minimum into the basin of attraction of another min- 
imum. The subsequent local geometry optimization part 
brings us then into the local minimum of this basin of 
attraction. In the original publication 6] it was already 
pointed out that the BEP principle is at least partly re- 
sponsible for the success of the minima hopping method. 
If the MD trajectory has a small kinetic energy Ekin, it 
can not go over very high barriers and it is thus more 
likely to reach the basins of attraction of low lying min- 
ima. It was shown that the number of local minima that 
was visited before the global minimum was found de- 
creases when the kinetic energy Ekin of the trajectory is 
reduced. Fig. [7J demonstrates the MDBEP principle for 
the Lennard-Jones cluster consisting of 55 atoms. There 
is a very strong correlation between the energy of the 
MD trajectory and the number of minima that are vis- 
ited before the global minimum is found. The data for 
Fig. [7J and the following similar figure were obtained by 
performing MHM runs that are stopped once the global 
minimum is found for different but fixed kinetic energies 
Ekin (he. 0i = 02 = 0a = 1 using the notation of ref. [3]) 
in a reasonably chosen energy interval. Subsequently we 
plot the values of Ekin versus the number of local minima 
that were visited before the global minimum was found. 
The potential energy of the local minimum from which 
the MD trajectory starts is set to zero. In this way the 
kinetic energy is the total energy of the MD trajectory 
and by energetic reasons it can not cross barriers higher 
than Ekin relative the starting minimum. Only new and 
accepted local minima are counted. In order to achieve 
better statistics we perform for each fixed Ekin 100 MHM 
runs (for Fig.[7jthe average is taken over 1000 runs), and 
take for the plots the averaged number of visited local 
minima. 

Since the Lennard Jones potential is a drastic simplifi- 
cation of the true inter-atomic interactions in solid state 
systems one might wonder whether the MDBEP principle 
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FIG. 7: The MDBEP principle for the Lennard- Jones cluster 
of 55 atoms. 
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FIG. 9: The MDBEP principle for the Morse cluster of 38 
atoms with p = 10.0 



also holds true for more realistic interactions. Using the 
MHM, we will therefore examine in the following the va- 
lidity of the MDBEP principle for other systems, namely 
Morse clusters and silicon clusters described both by a 
force field and a tight binding scheme. 

Fig. [5] and Fig. [5] present our results for the Morse clus- 
ters of 38 atoms with p — 6.0 and p = 10.0. Large values 
of p lead to a interaction that varies over shorter length 
scales. As a consequence the potential energy surface 
becomes more rugged and has significantly more local 
minima [3]. As a consequence considerably more min- 
ima are visited before the global minimum is found. The 
MDBEP principle is however well observed in both cases. 
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FIG. 8: The MDBEP principle for the Morse cluster cluster 
of 38 atoms with p = 6.0 

Fig. [TOland Fig.[TT]present our results for the Si™ clus- 
ter within the Lenosky tight binding scheme [lfj and 
for the S133 cluster within the Lenosky force field fill ] . 
In contrast to the Lennard Jones and Morse potentials 
the silicon force field has much more complicated inter- 



actions that depend not only on the distance between 
atoms but also on things like the bond angles. Tight 
binding schemes are the simplest way to treat solid state 
systems at a quantum mechanical level. The Lenosky 
tight binding scheme gave a very good agreement with 
the DFT energies |7| and can be considered as a reliable 
approximation to a precise density functional treatment 
of silicon clusters. The MDBEP principle is valid in both 
cases which demonstrates that the MDBEP principle is 
also valid for realistic interactions and in particular for 
quantum mechanical interactions. 
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FIG. 10: The MDBEP principle for the Lenosky tight binding 
cluster of 20 atoms. 



The fact that for small values of Ekin the global min- 
imum is found after having visited only a small number 
of local minima does not imply that the computational 
time in the MHM is continuously decreasing with smaller 
values of Ekin ■ If Ekin is getting too small the system has 
to make many attempts before succeeding to escape from 
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FIG. 11: The MDBEP principle for the Lenosky force field 
cluster of 33 atoms. 



the basin of attraction of the current minimum and this 
will actually lead to an increase in the computer time. 

In practice, the shortest computation time can be ob- 
tained by giving the MD trajectories initial velocities 
that have large components in the subspace of low cur- 
vature of the Hessian matrix. Due to the fact that low 
energy saddle points often lie at the end of low-curvature 
modes 17, EH one can in this way even with low 
energy trajectories rapidly escape from the present min- 
imum. 

In summary, we have shown that the BEP principle 
can be extended to MD trajectories. We call this ex- 
tended principle MDBEP principle and it says that MD 
trajectories with low energy are more likely to lead into 
basins of attraction of low energy configurations. Having 
verified this principle numerically for several systems we 
believe that it is valid for any solid state system. 

We thank the Swiss National Science Foundation for 
the financial support of our research. 
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